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1  INTRODUCTION 


Aircraft  structures  are  traditionally  passive  structures  optimized  for  a  given  flight  scenario  at 
the  minimal  cost  or  weight.  During  flight  these  structures  perform  efficiently  at  the  flight 
segment  it  was  optimized  for  but  at  less  optimal  efficiency  at  other  flight  segments.  This 
design  trait  limits  the  aerodynamic  efficiency  and  cost  effectiveness  in  multi-mission 
operations  and  obstructs  the  aircraft’s  versatility.  Ideally,  it  would  be  desirable  to  engineer 
aircraft  structures  that  can  operate  more  efficiently  across  the  entire  flight  envelope  with 
minimal  compromise.  This  inevitably  demands  transformation  in  the  size  and  shape  of 
aircraft  geometry  to  either  adapt  to  the  dynamic  aerodynamic  characteristics  or  to  induce 
favourable  aerodynamic  characteristics.  Therefore  based  on  the  airfoil  being  the 
aerodynamically  predominate  component  of  the  entire  aircraft,  there  is  the  need  for  a 
geometrically  adaptive  airfoil  (or  smart  wing  and  morphing  wing  as  cited  in  other  literatures) 
to  realize  the  improvement  in  efficiency.  It  has  always  been  the  general  consensus  of 
researchers  that  at  least  in  theory  an  adaptive  airfoil  is  more  effective  and  efficient  than  a 
passive  airfoil.  From  an  idealistic  aerodynamic  standpoint,  an  adaptive  airfoil  benefit  the 
aircraft  by  enhancing  L/D  ratio,  wake  control,  and  stall  characteristic  through  alterations  in 
its  airfoil  profile.  From  an  operational  viewpoint,  an  adaptive  airfoil  benefits  the  aviation 
industry  through  fuel  savings,  carbon  emission  reduction,  and  allows  greater  utilization  of 
aircraft  through  expansion  of  mission  profile.  One  such  analytic  study  (Bolonkin  and  Gilyard 
1999)  concluded  that  for  a  subsonic  transport  with  an  adaptive  airfoil  that  has  the  ability  to 
exert  trailing  edge  camber  deflection  can  significantly  reduce  drag  across  its  flight  envelope, 
especially  for  non-standard  flight  conditions  (up  to  10%)  while  less  for  cruise  (3%).  The 
reduction  in  drag  can  be  translated  into  monetary  term  through  fuel  savings  and  is  valued  at 
$US300,000  per  annum  for  a  transport  aircraft. 

The  design  of  adaptive  airfoil  remains  a  difficult  challenge  despite  recent  research  efforts  in 
both  industry  and  academic  arena.  This  can  possibly  be  attributed  to  two  compounding 
factors;  the  inherent  complexity  of  the  problem  and  the  forms  of  approach  taken  to  tackle 
them.  First  of  all,  the  design  objectives  in  adaptive  airfoil  are  competing  and  contradicting. 
Airfoils  are  traditionally  passive,  rigid  load  bearing  structures  and  by  infusing  compliance  or 
adaptiveness  into  the  design  environment  the  structure  will  inevitably  compromise  its 
stiffness,  thus  weaken  its  load  bearing  capability.  This  brings  the  difficult  task  of  striking  the 
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right  balance  between  stiffness  and  adaptiveness,  the  airfoil  needs  to  be  sufficiently  stiff  to 
sustain  the  aerodynamic  loads  and  structural  loads  without  suffering  un-wanted  shape 
change,  while  simultaneously  exhibits  features  of  flexibility  to  carry  out  controlled  shape 
adaptation.  Secondly,  traditional  airfoil  design  approach  has  been  largely  based  on  drawing 
experience  from  previous  designs,  and  this  design  approach  seems  to  have  resulted  in  similar 
stable  airfoil  structural  configurations,  as  evidenced  by  the  wing  design  in  the  last  50  years 
of  passenger  jet  aircraft.  This  “rely  on  experience”  design  approach  has  sufficed  until  now 
for  the  passive  airfoil  designs,  but  with  the  newly  added  need  for  adaptiveness  mentioned 
previously  the  design  environment  has  changed  drastically.  The  conventional  arrangement  of 
spars,  ribs,  and  stringer  will  not  suffice  the  design  needs,  and  with  limited  history  in  adaptive 
airfoil  operation  little  experience  can  be  drawn  when  designing  one.  This  has  led  to  many 
design  proposals  derived  from  ad-hoc  approaches,  and  these  approaches  are  mostly  resulted 
from  trial  and  error.  The  design  objectives  are  usually  qualitative  and  ambiguous,  which  in 
turn  causes  difficulty  when  comparing  the  performance  against  the  original  goal  or  against 
another  design,  plus  further  improvements  to  the  design  most  likely  to  be  decided  in  ad-hoc 
fashion  also. 

To  realize  the  adaptive  airfoil  design,  an  alternative,  more  systematic,  and  more  theoretical 
design  approach  is  desired.  This  report  brings  a  new  design  philosophy  called  “unit  cell 
approach”  to  envision  new  adaptive  airfoil  configuration,  furthermore  topology  optimization 
was  chosen  as  the  computing  tool  to  generate  the  design  in  a  systematic  procedure  under 
computer  aided  design  environment.  Although  this  research  revolves  around  the  issue  of 
adaptive  airfoil  the  implications  of  this  work  are  much  broader,  any  vehicle  that  operates 
across  multiple  environment  settings  and  in  different  operation  modes  can  potentially  benefit 
from  infusing  shape  adaptation  ability  to  the  structure.  Therefore,  this  research  is  extremely 
important  to  the  design  of  current  and  future  generations  of  vehicles  alike. 

The  specific  adaptation  feature  in  study  is  adaptive  camber  intended  at  the  leading  edge  or 
trailing  edge  of  the  airfoil  for  maximum  aerodynamic  effect.  The  main  purpose  of  implanting 
adaptive  camber  at  the  leading  and  trailing  edge  is  to  derive  aerodynamic  efficiency  at  non¬ 
standard  (cruise)  flight  segments.  Numerical  examples  as  well  as  prototype  testing  will  be 
presented  to  demonstrate  a  range  of  issues  such  as  the  necessity  of  multi-objective 
considerations  and  high  fidelity  shape  control.  A  new  method  of  handling  multi-objective 
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topology  optimization  by  combining  SIMP  method  (Solid  Iso-tropic  Material  with 
Penalization)  with  Physical  Programming  (PP)  is  proposed  and  demonstrated  in  this  research. 

Secondly,  alternative  topology  optimization  method  will  he  investigated  in  order  to  reduce 
the  computation  cost.  A  novel  topology  optimization  method  termed  “Evolutionary  level  set” 
(ELS)  is  developed  and  demonstrated  with  numerical  examples.  The  method  of  ELS  is 
unique  for  taking  advantage  of  the  computational  effective  bio-evolutionary  principles  of 
ESO  and  also  the  implicit  free  boundary  representation  of  LSM.  The  phase  interface  between 
structural  boundary  and  void  is  numerically  described  at  each  iteration  in  implicit  fashion  as 
dynamic  level  sets  of  a  sequence  of  iso-surfaces.  An  evolutionary  algorithm  is  developed  to 
advance  the  movement  of  the  phase  interface  and  update  the  level  set  function.  As  a  result 
the  solving  of  the  Hamilton-Jacobi  PDEs  through  explicit  time-marching  schemes  of 
conventional  LSM,  together  with  the  complex  shape  derivative  analyses  as  well  as  the 
numerical  difficulties  are  bypassed.  Simultaneously  the  major  advantages  of  the  level  set 
based  boundary  representation  seheme  are  maintained. 


2  DESIGNING  ADAPTIVE  AIRFOIL  VIA  SIMP-PP 

2. 1  Unit  cell  philosophy 

A  unit  cell  structure  is  a  building  block  used  in  the  assembly  of  a  larger,  global  structure 
through  repeatedly  linked  networks.  It  is  required  to  carry  the  generic  performance  feature  of 
the  global  strueture  (sueh  as  shape  adaptation)  and  by  skilled  design  this  feature  can  be 
accumulated  and  magnified  over  a  network  of  repeatedly  linked  cells.  Certain  unit  cell 
struetures  sueh  as  the  Kagome  truss  (Symons,  Hutehinson  et  al.  2005)  holds  potential  in  the 
design  of  morphing  structures  that  require  the  structure  to  behave  in  isotropic  manner  under 
external  loading  while  possess  the  ability  to  carry  out  large  deformation  under  internal 
actuation.  When  it  comes  to  shape  adaptive  structures  the  implication  of  a  unit  cell  network 
design  have  several  benefits,  first  of  all  the  adaptation  is  generated  through  the  accumulation 
of  each  individual  adaptive  cell,  meaning  the  adaptation  is  distributed  across  large  regions  of 
the  structure  and  not  concentrated  at  a  local  region,  this  implies  that  the  individual  shape 
change  of  each  cell  does  not  neeessarily  have  to  be  forceful,  thus  lowering  both  the  power 
demand  of  the  actuator  which  possibly  carries  a  positive  impaet  on  the  size  of  the  actuator 
and  also  the  stress  level  on  the  strueture  which  widens  the  material  selection  options  for  the 
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host  structure.  Secondly,  by  distributing  the  adaptive  control  across  a  large  region  while  each 
unit  cell  can  be  actuated  independently,  this  effectively  makes  each  unit  cell  an  independent 
control  surface;  the  possible  combination  of  actuation  control  setting  thus  increases 
exponentially  with  the  number  of  unit  cells,  therefore  allowing  the  adaptive  network  to 
achieve  complex,  non-linear  profiles  with  high  fidelity. 

This  research  aims  to  present  a  systematic  design  method  for  the  adaptive  camber  airfoil 
structure  inclusive  of  the  actuation  material  by  using  the  method  of  topology  optimization  to 
design  an  adaptive  unit  cell  which  has  the  adaptation  features  representative  of  camber 
deflection,  and  then  apply  the  concept  of  unit  cell  network  to  construct  the  sizable  airfoil 
structure  from  the  assemblage  of  multiple  unit  cells.  This  approach  allows  the  simultaneous 
design  of  the  airfoil  structural  geometry  as  well  as  the  size,  shape,  form  and  distribution  of 
the  actuators  for  the  purposes  of  shape  adaptation  and  load  bearing. 

In  this  research,  a  simple  unit  cell  network  is  proposed  as  a  starting  point,  the  unit  cell 
network  is  created  by  repeatedly  linking  unit  cells  in  series.  For  the  purpose  of  camber 
adaptation  at  the  trailing  edge,  the  unit  cell  network  is  attached  to  the  rear  spar,  and  is 
tapered  to  accommodate  the  airfoil  geometry  (shown  in  Figure  1).  The  rear  spar  provides  the 
rigid  support  to  nullify  the  reaction  forces  caused  by  the  actuation. 


Figure  1  Unit  cell  network  illustration 
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The  objectives  in  the  design  of  the  unit  cell  structure  involve  considering  both  structural 
stiffness  and  adaptation.  In  the  case  of  structural  stiffness  the  unit  cell  is  designed  to  carry 
(for  a  2-D  case)  transverse  loading  transmitted  from  the  wing  skin  and  bending  moments 
along  the  chord- wise  direction  (see  Figure  2).  As  for  adaptation  the  unit  cell  is  designed  to 
achieve  positive  or  negative  camber  deflection  by  repositioning  the  top  and  bottom  comers 
of  its  rear  edge  in  opposite  directions,  thus  deflecting  the  remaining  sections  of  the  unit  cell 
network.  Each  unit  cell  in  the  network  is  controlled  independently,  which  allows  for  the 
accumulation  of  camber  deflection  over  the  consecutive  cells,  leading  to  a  substantial 
camber  deflection  reaching  the  wing  tip.  However,  the  adaptive  camber  stmcture  can  also 
take  non-monotonous  camber  deflection  by  simply  reversing  the  actuation  motion  at  a 
certain  unit  cells  inside  the  network  to  enable  the  camber  acquiring  a  profile  of  a  high-order 
polynomial. 


2.2  Multi-objective  topology  optimization  method  SIMP-PP 

As  mentioned  previously  one  of  the  challenging  aspects  in  adaptive  airfoil  design  is 
achieving  a  delicate  balance  between  design  objectives;  namely  stiffness  and  flexibility. 
Thus  the  utilization  of  multi-objective  optimization  scheme  is  a  necessity  for  a  reliable 
stmcture  that  can  perform  adequate  adaptation  while  resist  aerodynamic  loading.  The 
advantages  of  SIMP-PP  compared  to  SIMP  is  demonstrated  in  the  resulting  topology 
optimization  comparisons  later.  Due  to  the  fact  the  theory  of  SIMP-PP  was  thoroughly  laid 
out  in  the  previous  annual  report,  it  shall  not  be  further  explained  in  this  section.  For  more 
information  on  SIMP-PP  the  reader  is  referred  to  the  previous  report  or  the  paper  (Lin,  Luo 
et  al.  2009). 
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2.3  Multi-material  interpolation  scheme 


In  an  arguably  ad-hoc  fashion,  a  logical  but  empirical  design  parameter  was  made.  The 
adaptive  structure  is  designed  to  constitute  multi-material  structure.  With  at  least  one  type  of 
passive  material  acting  as  the  host  structure,  and  one  type  of  active  material  acting  as  the 
actuation  material.  In  this  particular  study  a  bi-material  compound  is  designed  with  topology 
optimization.  This  study  employs  a  hybrid  multi-material  interpolation  scheme  (Bendsoe  and 
Sigmund  1999;  Sigmund  2001)  which  combines  the  standard  power  law  method  with  the 
original  Hashin-Shtrikman  bound  interpolation  (Hashin  and  Shtrikman  1963).  Trials  have 
shown  that  bi-material  designs  offer  better  performance  than  the  single  material  design;  this 
is  mostly  due  to  the  fact  that  bi-material  design  takes  advantage  of  the  different  material 
properties  and  enables  unimorph  configuration  which  is  naturally  suited  for  certain  structural 
adaptation,  i.e.  deflection.  Single  material  compliant  mechanisms  on  the  other  hand  typically 
rely  on  pure  tension-compression  mechanisms  to  derive  the  desired  adaptation  and  can  be 
less  effective  in  some  cases.  There  are  effectively  three  distinct  material  phases  in  this  bi¬ 
material  design,  two  solid  types  and  void.  It  is  necessary  to  interpolate  between  the  phases 
during  the  optimization  so  that  the  effective  material  properties  such  as  elastic  modulus  E, 
thermal  expansion  coefficient  a,  electrical  and  thermal  conductivity  coefficient  |3  and  j  are 
used  in  the  solving  of  governing  equations.  With  this  approach  the  material  properties  are 
interpolated  as 

^(x,  v)  =  v"  [(1  -  /c)<^f  (x)  +  (x)] ,  0  <  /c  <  1 

g?(x,  v)  =  [(1  -  a:)^®  (x)  -r  /c(p^^  (x)] ,  0  <  x  <  1 

/?(x,  V)  =  V"  [(1  -  ^)/?f  (x)  +  (X)] ,  0  <  V  <  1 

Z(x,  V)  =  V"  [(1  -  W  +  A-jf  (x)] 

Parameters  such  as  and  are  the  upper  and  lower  H-S  bound  value  for  the  related 

parameter,  k  is  the  penalization  factor  between  the  upper  and  lower  bound  values.  ^  and  (p 


are  the  bulk  and  shear  modulus  which  relates  to  the  effective  elastic  modulus  E  by  Eq  (2)  for 
a  2-D  problem,  x  and  v  are  the  two  variables  for  which  the  material  properties  of  each 
elements  depends  upon. 


K{x,v)  +  G(x,v) 


(2) 
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According  to  the  theory  by  Hashin-Shtrikman  (1963),  the  thermal  expansion  coefficient  is 
slightly  different  as  it  relies  only  on  the  bulk  modulus  and  is  given  as 


^  _  ^1^1  (^2  ^l)  ^2^2  (^1  ^2) 


(3) 


where  parameters  with  subscripts  1  and  2  are  the  intrinsic  material  properties  for  solid 
material  1  and  2  respectively.  This  report  only  presents  the  interpolation  scheme  for  the 
electro-thermal-mechanical  related  material  properties,  and  the  details  for  resolving  the 
various  upper  and  lower  H-S  bound  can  be  found  in  (Hashin  and  Shtrikman  1963;  Sigmund 
2001). 


2.4  Adaptive  airfoil  Formulation 


Writing  the  problem  formulation  for  topology  optimization  in  the  Karush-Kuhn-Tucker  form, 
we  have 

Minimize :  f(x) 

'  [g.(x)>0  j  =  (4) 

Sub  ect  to:  \ 

[h.(x)  =  0  i  =  l,...,n 


with  f(x)  being  the  objeetive  funetion,  g/x)  the  inequality  constraint,  and  hi(x)  the  equality 
constraint. 


The  adaptive  airfoil  formulation  can  be  stated  as 

'  Minimize^ :  fSX)  =  g(e{X),J {X)) 

\Vx,V2,...,Vn  ) 


Subject  to: 


Vjc  V  -K  <0, 

/  j  e  e  1  ’ 

e=l 

X(x^-i)v^-r,<o, 

e=l 

v-1  G  Q  =  0, 

<  X  G  Q  =  0, 

V  G  Q  =  0, 


0<x^<l,0<v^<l, 
Equalibrium  equations 


(5) 
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fa(X)  is  the  aggregate  function  of  the  multiply  objective  functions  calculated  based  on  the 
deflection  angle  0(X)  and  the  strain  energy  J(X).  The  deflection  angle  indicates  the 
adaptation  performance  of  the  unit  cell  while  the  strain  energy  serves  as  the  classic  measure 
of  stiffness.  The  aggregate  function  fa(x)  is  calculated  by  the  method  of  physical 
programming.  The  design  variable  vector  X  consists  of  variables  and  V],...,vn,  with 

N  being  the  number  of  elements  meshing  the  design  domain.  Variable  x  interprets  between 
the  passive  material  and  the  solid  material,  variable  v  interprets  between  solid  material  and 
void.  A  clear  illustration  is  given  in  Table  1  below. 

Table  1  Multi-material  interpolation  method 


V=1 

v=0 

X=1 

Solid  region  with  active  material 

Void  region 

x=0 

Solid  region  with  passive  material 

Void  region 

The  parameters  Vj  and  V2  in  the  two  inequality  eonstraints  are  the  maximum  volume 
eonstraints  for  the  aetive  and  passive  materials  respeetively.  In  most  topology  optimization 
routines  the  design  domains  are  generally  of  reetangular  dimensions.  To  enforee  a  tapered 
airfoil  strueture  fixed  solid  and  void  regions  are  ineorporated  into  the  design  domain.  The 
parameter  Q  in  the  equality  eonstraint  equations  represents  the  passive  region  inside  the 
design  domain  as  was  shown  in  Figure  3.  The  parameter  Q  represents  the  void  region 
neeessary  to  enforee  a  tapered  design  spaee  inside  a  reetangular  domain. 
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The  equilibrium  equations  refer  to  the  governing  equations  of  the  system.  The  aetuation  of 
most  adaptive  eore  is  through  eertain  external  stimuli  whieh  the  aetuator  material  is  sensitive 
to,  and  respond  by  indueing  meehanieal  energy  upon  the  passive  host  strueture.  The 
orthodox  aetuator  in  adaptive  eore  design  is  PZT  material,  whieh  is  based  on  an  eleetro- 
meehanieal  system.  The  aetuation  system  defined  in  this  work  is  an  example  of  eleetro- 
thermo-meehanieal  system  indueed  by  uniform  Joule  heating,  although  any  other  aetuation 
method  ean  be  formulated  in  similar  manner.  Multiphysies  aetuation  systems  sueh  as  this 
one  have  been  widely  studied  reeently  due  to  many  of  its  advantages,  whieh  inelude 
relatively  large  shape  adaptation  for  given  aetuation  energy,  good  eontrollability,  fabrieation 
ease  and  amenability  (Sigmund  2001;  Luo,  Tong  et  al.  2009).  A  general  model  for  eleetro- 
thermo-meehanieally  loaded  struetures  is  to  sequentially  solve  a  set  of  eoupled  boundary 
value  problems  that  span  the  eleetrieal,  thermal  and  elastie  domains  as  shown  in  Eq  (6).  By 
whieh  an  eleetrieal  input  is  indueed  upon  the  system  to  trigger  a  thermal  response  (typieal  of 
Joule  heating),  and  the  thermal  response  eonsequently  results  in  a  meehanieal  output  due  to 
the  intrinsie  properties  of  the  host  strueture  whieh  serves  as  an  aetuation  foree  to  drive  the 
deformation.  This  system  is  eonsidered  to  be  weakly  eoupled,  whieh  denotes  that  the 
eleetrieal  equation  is  independent  of  the  heat  equation,  and  the  heat  equation  is  unattaehed  to 
the  elastie  equation  in  a  similar  way. 

Eleetrieal  equation:  (x,  v)U^  (v,  v)  = 

Thermal  equation:  (v,  v)U^^  (v,  v)  =  (U^  (v,  v),  v,  v)  (6) 

Meehanieal  equation:  (x,  (v,  v)  =  F^^^  (v,  v),  v,  v) 

The  sub  index  /,  II  and  III  represent  the  eleetrieal  system,  thermal  system  and  meehanieal 
system  respeetively.  K/,  U/,  and  F/  are  the  global  eleetrieal  eonduetivity  matrix,  voltage  field 
and  eleetrieal  load  veetor.  K//,  U//,  and  F//  are  the  thermal  eonduetivity  matrix,  temperature 
field  and  thermal  load  veetor.  K///,  U///,  and  F///  are  the  global  stiffness  matrix,  displaeement 
field  and  loading  veetor.  The  three  sets  of  equations  will  need  to  be  solved  in  sequenee 
starting  with  eleetrieal  equation,  then  thermal  and  finally  meehanieal  to  represent  the 
idealized  ease  of  eleetrieally  driven  heating  and  then  the  meehanieal  effeet  of  the  thermal 
strain. 
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2.5  Case  study 


This  section  presents  the  case  study  which  lead  to  the  prototype  design  in  chapter  3,  The 
design  domain  in  this  case  is  a  rectangular  domain  100x80  mm  with  1mm  thickness,  meshed 
using  8000  quad4  plane  elements.  The  left  edge  of  the  unit  cell  is  fixed  and  a  pre-defined 
solid  region  is  assigned  to  the  right  edge  to  represent  the  cell  wall  which  is  to  he  deflected. 
The  taper  ratio  is  set  to  0.9.  The  loading  condition  on  the  unit  cell  are  the  transverse  loading 
transferred  to  the  unit  cell  structure  through  the  skin,  and  the  accompanying  bending  stress 
transmitted  through  the  adjacent  cell.  They  are  modelled  as  an  evenly  distributed  vertical 
force  and  a  linearly  distributed  horizontal  force  respectively.  The  displacement  output  ports 
Uouti  and  Uout2  are  at  the  top  and  bottom  comers  of  the  left  edge  and  in  opposite  directions  to 
impose  deflection. 

The  material  properties  of  this  bi-material  stmcture  is  set  to  have  one  material  being 
particularly  sensitive  to  electric  stimulation  with  a  small  volume  constraint  and  have  the 
other  material  being  inert  with  a  large  volume  constraint.  The  later  in  real  life  would  serve  as 
the  passive  host  stmcture  while  the  former  serve  as  the  embedded  active  material  or  actuator. 
It  is  assumed  (for  practical  reasons)  that  in  practical  applications  the  electrical  actuation  is 
applied  strictly  to  the  active  material  and  not  the  entire  stmcture  and  in  forms  of  either  a 
defined  voltage  or  current  from  a  power  supply.  In  this  study  the  loading  condition  is 
implemented  as  the  entire  domain  is  subjected  to  an  electrical  actuation  in  the  form  of  a 
uniform  voltage,  rather  than  at  specific  terminal  locations.  Even  though  in  this  loading 
condition  all  material  in  the  domain  will  be  under  electrical  actuation,  due  to  the  inert 
characteristics  of  the  passive  material  the  load  will  have  little  or  no  effect  on  the  passive 
material,  hence  the  underlining  design  concept  of  the  actuation  source  comes  only  from  the 
active  material  is  still  valid,  and  only  the  actuator  needs  to  be  powered  in  actual  application. 

The  passive  host  material  is  titanium  (E  =  150GPa,  v  =  0.31,  a  =  8.6e-6K-l,  p  =  2.38e6Q- 
lm-1,  X  =  21.9Wm-lK-l)  and  the  active  material  is  an  artificial  actuator  with  a  high  thermal 
expansion  coefficient  (E  =  15GPa,  v  =  0.31,  a  =  -2e-3K-l,  P  =  3e6Q-lm-l,  %  =  50Wm-lK- 
1)  with  volume  fractions  of  15%  and  3%  of  the  total  design  domain  respectively.  Topology 
optimization  was  conducted  multiple  times  by  tweaking  the  relative  importance  of  the  two 
objectives  to  observe  the  difference  between  a  stiffness  intensive  optimization  and  a 
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performance  intensive  optimization.  For  the  sake  of  numerical  simplicity  linear  elastic 
structures  are  analysed. 

The  articulations  of  preferences  are  done  in  the  same  fashion  as  outlined  previously  (Lin, 
Luo  et  al  2009),  and  in  this  case  the  preferences  for  the  stiffness  are  held  constant  in  each 
routine  while  the  adaptation  preferences  are  altered.  The  articulation  for  the  adaptation 
preferences  is  displayed  in  Figure  4.  The  optimization  results  for  the  compliant  unit  cell  are 
presented  in  the  following  section;  Table  2  includes  the  optimal  topologies,  the  stiffness 
intensive  to  adaptation  intensive  optimizations  are  arranged  from  case  a  to  f,  with  case  a 
being  the  most  stiffness  intensive  for  carrying  the  defined  loads  and  case  f  the  most 
adaptation  intensive.  The  grey  areas  represent  the  active  material  while  the  black  areas 
represent  the  passive  material.  Table  3  and  Table  4  shows  the  objective  functions  summary. 


Figure  4  Articulation  of  adaptation  preferences  in  cases  a)  to  f) 
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Table  3  SIMP-PP  optimization  results  for  tapered  unit  cell 


Case  Operation  Highly  Des.  Desirable  Tolerable  Undesirable  Highly  Undes. 


fa 

fa 

fa 

fi4 

fi5 

a 

maximize  camber 
(deg) 

4.0 

1.0 

0.6 

0.45  ^ 

0.4 

minimize  SE 

(j) 

2.5 

- 

12.5 

32.5 

55.5 

b 

maximize  camber 
(deg) 

4.0 

2.2 

- 

0.6 

0.4 

minimize  SE 

(j) 

2.5 

7.5 

p.5 

22.5 

55.5 

Q 

maximize  camber 
(deg) 

4.0 

3.0 

1 

1 

1 

1 

1 

1 

.....0-4 

minimize  SE 

(J) 

2.5 

7.5 

- 

22.5 

55.5 

d 

maximize  camber 
(deg) 

4.0 

3.2 

1.6 

0.4 

minimize  SE 

(J) 

2.5 

7.5 

- 12A^ 

22.5 

55.5 

maximize  camber 

4.0 

3.5 

2.3 

0.4 

e 

(deg) 

minimize  SE 

2.5 

7.5 

^ - 

12.5 

^22.5 

55.5 

(J) 

► 

maximize  camber 

4.0 

3.9 

3.0 

0.3 

f 

(deg) 

minimize  SE 

2.5 

7.5 

12.5 

22J 

55.5 

(Ji 
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Table  4  Objective  function  summary  of  tapered  unit  cell 

Case_ SE  (J) Output  (deg) 


a 

8 

0.4 

b 

11 

1.6 

c 

14 

2.2 

d 

18 

2.4 

e 

26 

3.1 

f 

37 

3.7 

From  both  the  performance  summary  and  the  topological  results,  it  can  be  deduced  that 
when  the  articulation  of  preferences  favours  stiffness  over  adaptation,  the  topology  is 
characterised  by  structural  features  such  as  hinge-less,  thick  structural  members,  truss 
networks,  and  axle  actuation  with  either  tension  or  compression  mechanisms.  Whereas  when 
the  articulation  of  preferences  favours  adaptation  over  stiffness,  the  topology  is  often 
characterised  by  features  such  as  hinges  and  single  element  connections,  complex  uni-morph 
mechanisms,  and  complicated  structural  arrangement  with  large  number  of  structural 
members.  In  fact,  the  same  observations  were  made  previously  in  the  design  of  rectangular 
adaptive  core,  which  did  not  involve  multi-objective  optimization  with  stiffness  as  a  design 
objective. 


2.6  Summary 

The  research  presented  in  this  section  illustrates  the  importance  of  a  well  formulated 
optimization  routine  to  the  design  of  adaptive  airfoil  structure.  The  multi-objective  nature  of 
an  adaptive  airfoil  must  be  taken  into  account  in  order  to  produce  a  practical  and  applicable 
topology.  From  a  design  perspective  cases  a)  b)  and  c)  (which  all  has  preference  articulation 
favouring  stiffness)  are  much  more  desirable  than  others  due  to  their  structural  simplicity 
and  stability.  The  parameters  employed  in  this  report  serve  as  a  guide  for  future  studies  and 
does  not  need  to  be  followed  to  the  letter.  Indeed  it  would  be  very  difficult  to  find  matching 
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actuator  material  that  precisely  matches  the  one  defined  in  the  optimization,  or  matching 
actuator  with  the  exact  size  and  functionality.  Rather  topology  optimization  serves  as  a  good 
starting  point  to  explore  different  design  options,  and  it  should  be  up  to  the  designer  to 
further  enhance  the  topology  in  the  later  stages  of  the  design  process  to  consider  for  other 
design  objectives  such  as  manufacturing,  cost  and  maintenance. 

Comparing  the  design  selected  for  prototyping  from  this  report  with  the  previous  one  (see 
Figure  5).  The  current  design  is  by  far  possibly  the  most  realistic  design  option.  It  has  several 
key  benefits;  firstly  its  structure  distribution  is  balanced  (for  the  passive  material),  one  could 
easily  enhance  the  actuation  force  of  the  system  by  adding  an  additional  actuator  at  the  lower 
half  of  the  structure.  Secondly,  its  structure  is  relatively  easy  to  manufacture  due  to  its 
simple  and  novel  layout.  In  fact,  the  second  prototype  was  manufactured  off  the  topology 
with  little  to  no  post  processing,  where  as  for  the  previous  one  post  processing  was  necessary 
to  reform  the  single  point  hinges  in  the  system  to  compliant  hinges.  Thirdly,  the  actuation 
system  depicted  in  the  topology  is  easy  to  implement,  the  topology  depicted  a  single  linear 
actuator  which  is  common  in  industry.  One  notable  difficulty  encountered  in  this  research 
was  to  find  a  suitable  actuator  that  was  appropriate  for  the  size  of  the  prototype. 


Figure  5  Prototyped  unit  cells;  left)  in  current  report,  right)  previous  report 


3  TESTING 

3. 1  Design  selection 

The  topology  employed  for  prototyping  and  testing  was  chosen  as  case  b)  from  Table  2.  The 
topology  was  analysed  in  FEA  to  determine  its  structural  performance  and  material 
requirements  in  both  single  cell  mode,  and  three  cell  network  mode.  In  summary  the 
selection  process  was  conducted  as  follows 
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1 .  Select  optimal  topology 

2.  Selection  manufacturing  method  and  passive  material 

3.  Analyse  actuation  energy  required  based  on  the  passive  material  through  FEA 

4.  Select  actuator  based  on  the  actuation  energy 

5.  Adjust  the  design  (topology)  if  necessary  to  accommodate  the  actuator 

Due  to  the  size  and  geometry  of  the  intended  prototype  and  its  adaptation  feature,  the  most 
realistic  manufacturing  method  is  rapid  prototyping  using  flexible  polycarbonate  material. 
Rapid  prototyping  would  provide  a  reasonably  accurate  complex  3-D  product  at  a  short  time 
frame  with  minimal  material  wastage.  The  material  used  in  rapid  prototyping  is  industrial 
polycarbonate,  a  common  component  among  automotive  and  aerospace  design.  It  has  an 
elastic  modulus  of  2GPa,  tensile  strength  of  52]V[Pa  and  a  yield  strain  of  3%.  The  exhibition 
of  a  high  yield  strain  and  superior  mechanical  properties  compared  to  other  thermoplastics 
makes  a  strong  case  for  the  employment  of  polycarbonate  in  this  lab  testing. 

FEA  tests  were  conducted  to  analyse  the  actuation  energy  requirement  on  the  host  structure 
made  of  polycarbonate  and  the  resulting  adaptation  performance.  The  analysis  was  done  on  a 
single  tapered  unit  cell  and  a  unit  cell  network  made  of  3  cells.  The  analysis  is  displayed  on 
the  figures  and  graphs  below.  Figure  6  and  Figure  7  display  the  deformed  unit  cell/network 
with  Von-Mises  strain  contour.  In  the  case  of  Figure  7  it  is  assumed  that  the  actuators  can 
provide  sufficient  blocking  force  to  not  only  hold  the  adaptive  structure  at  a  given  shape 
configuration,  but  also  prevent  the  actuation  of  neighbouring  cells  to  interfere  with  the 
adaptation  of  the  current  one.  This  assumption  is  under  ideal  circumstances;  as  it  would  be 
difficult  to  estimate  the  degree  of  interference  of  neighbouring  actuation  have  upon  the 
current  cell. 
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Figure  7  Undeformed  and  deformed  network  with  strain  contour 


Figure  8  illustrates  the  relation  between  adaptation  achieved  and  the  actuation  energy 
required  for  the  first  unit  cell  under  adaptation  with  the  dimensions  150x150x5  and  a  taper 
ratio  of  0.9.  While  Figure  9  displays  the  relationship  between  strain  and  actuation  energy. 
Analysis  indicates  that  elevation  in  actuation  energy  increases  the  camber  deflection  with 
almost  linear  proportionality,  but  simultaneously  increases  the  strain  in  the  substructure. 
Which  may  lead  to  yielding  given  the  substructure  is  aluminium  or  typical  aerospace  grade 
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alloy.  While  it  is  not  expected  for  a  single  unit-cell  of  substructure  to  carry  out  camber 
deflection  larger  than  5  degrees  (flow  separation  will  start  to  occur  at  moderate  degrees  of 
camber),  the  material  selection  against  yielding  may  well  be  an  important  factor  in  any 
practical  application. 


Figure  8  Camber  deflection  vs  actuation  energy  for  single  cell  adaptation 


Figure  9  Von-Mises  strain  vs  actuation  energy  for  single  cell  adaptation 

The  aim  for  the  prototype  is  to  generate  a  total  trailing  edge  deflection  of  10  degree  within  a 
three  cell  tapered  network  as  indicated  in  Figure  10.  Hypothetically  if  each  cell  achieves  an 
average  deflection  of  three  degrees,  that  would  put  the  strain  at  0.7%,  which  is  well  below 
the  yield  strain  limit  of  the  polycarbonate  material. 


It  is  decided  to  employ  servo  motors  as  the  actuator  in  this  prototype  test  rather  than  using 
Nitinol  spring  actuators  in  the  previous  report.  The  servo  motors  -  although  heavier  than  the 
Nitinol  actuators,  have  several  attractive  attributes.  Firstly  servo  motors  have  excellent 
respond  time  given  they  are  completely  electrically  actuated,  whereas  in  the  case  of  Nitinol 
actuators  the  actuation  suffers  a  lag  period  for  the  Joule’s  heating  to  take  place.  Secondly  the 
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servo  motors  provide  exeellent  bloeking  foree  that  loeks  in  the  servo  in  the  aetuated  position 
onee  the  aetuation  load  is  removed.  And  thirdly  the  servo  motor  arm  has  suffieient  stiffness 
that  ean  truly  eontribute  to  the  load  bearing  eapability  of  the  adaptive  system  against  external 
loading,  making  this  an  integrated  aetuation  system  that  serves  not  only  as  a  souree  of 
aetuation  but  also  a  key  struetural  element.  The  downside  of  using  servo  motors  are  the 
inerease  in  aetuation  weight  (as  mentioned  above)  and  also  the  near  linear  but  not  strietly 
linear  motion  the  servo  provides.  The  weight  inerease  for  this  partieular  ease  is  eonsidered 
negligible  as  eaeh  servo  weighs  no  more  than  30g.  The  problem  with  the  non-linear  motion 
from  the  servo  was  minimized  by  linking  the  servo  to  the  aetuation  point  with  a  suffieiently 
long  arm  would  insure  the  rotational  motion  transfers  to  almost  linear  motion  at  the  other 
end  of  the  arm. 


3.2  3-D  projection 

Before  manufaeturing  takes  plaee  the  2-D  topology  was  partially  extruded  in  the  thiekness 
direetion  to  realize  a  3-D  strueture,  this  spanwise  extrusion  is  eommon  praetiee  in  adaptive 
airfoil  design.  The  extrusion  aehieves  a  final  partially  enelosed  thin  wall  strueture  that  helps 
to  stabilize  lateral  torsion,  and  reinforees  the  eelTs  edge  wall  to  resist  deformation  eaused  by 
the  aetuation.  The  strueture  was  further  modified  before  manufaeturing  to  aeeommodate  the 
aetuator  attaehments  and  rig  attaehments.  The  final  design  of  the  passive  eomponent  of  the 
adaptive  system  is  shown  in  Figure  11.  The  aetual  prototype  with  aetuators  attaehed  is 
shown  in  Figure  12,  eompared  to  the  original  topology  in  Table  2,  the  aetuation  system 
remains  almost  eompletely  similar  to  the  suggested  optimal  topology.  The  final  unit  eell 
network  has  the  dimension  379x66x180  mm^. 
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Figure  11  CAD  model  of  the  final  passive  component  in  the  adaptive  system 


Figure  12  Prototyped  adaptive  network 


3.3  Testing 

Testing  is  conducted  to  prove  and  showcase  the  advantages  of  employing  unit  cell  design. 

The  specific  aims  of  the  tests  are 

•  Illustrate  the  ability  of  the  unit  cell  design  to  accumulate  individual  cell  adaptation  to 
build  up  a  sizable  camber  deflection  of  10  degrees  over  the  three  cell  network 

•  Illustrate  that  the  adaptive  system  possesses  distributed  actuation  where  each  cell  can 
be  actuated  independently,  thus  making  each  cell  an  independent  control  surface 

•  Illustrate  the  high-fidelity  shape  adaptation  aspect  of  the  design  by  showcasing  the 
wide  range  of  airfoil  profiles  the  3  cell  unit  cell  system  can  bring 
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These  desired  aims  are  intended  to  highlight  the  benefit  of  the  design  approaeh/eoneept 
rather  than  highlight  the  aetual  qualitative  performanee  of  this  partieular  adaptive  system, 
sinee  the  qualitative  performanee  ean  easily  be  altered  by  using  different  material,  different 
aetuators,  and  different  aetuation  inputs.  The  tests  are  eondueted  by  plaeing  AP-1297 
position  sensor  probes  around  the  strueture  to  log  displaeement  ehange  at  speeifie  points  e.g. 
in  Figure  13;  the  logged  displaeements  ean  then  be  turned  into  defleetion  angles  by  referring 
to  the  referenee  position. 


Figure  13  Lab  setup 

In  the  first  test  ease  eondueted,  all  three  eells  were  aetuated  simultaneously  and  the  total 
defleetion  generated  was  traeked.  The  displaeement  plot  is  shown  in  Figure  15  and  the 
defleetion  plot  is  shown  in  Figure  16.  The  data  indieates  that  a  total  eamber  defleetion  of  9.4 
deg  was  reaehed  by  the  simultaneous  aetuation  of  all  three  eells  at  a  defleetion  rate  of  4.7 
deg  per  seeond.  The  aetual  defleetion  rate  was  mueh  faster  than  the  data  indieated  but  is 
disguised  by  the  noise/vibration  filtering  proeess.  In  both  plots  there  are  minor  vibration 
motions  during  the  “holding”  stage  and  after  the  retraeting  stage.  The  vibration  motion 
during  the  “holding”  stage  is  due  to  the  servo’s  bloeking  foree  at  holding  the  eamber 
defleetion  in  plaee.  The  vibration  motion  after  the  retraeting  stage  is  due  to  the  natural 
damping  in  the  host  polyearbonate  aeting  as  transient  residual  load  on  the  servos.  Figure  14 
shows  the  adaptive  eamber  strueture  before  and  after  aetuation. 
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Figure  14  Adaptive  camber  before  and  after  actuation 


Figure  15  Displacement  plot  of  simultaneous  actuation  of  adaptive  network 
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Time  (s) 


Figure  16  Deflection  plot  of  simultaneous  actuation  of  adaptive  network 


In  the  second  test  conducted,  the  three  cells  were  actuated  in  sequence  starting  from  the  first 
cell,  then  the  second,  then  the  third  under  the  same  actuation  input  to  achieve  deflection.  And 
once  after  the  total  deflection  has  been  reached,  the  unit  cells  are  de-activated  to  revert  back 
to  the  neutral  position  in  the  opposite  sequence  starting  with  the  third  cell,  then  the  second 
and  finally  the  first.  The  results  are  displayed  on  Table  6.  The  data  indicates  that  the 
adaptation  control  in  each  unit  cell  can  be  considered  as  independent  given  there  is  sufficient 
blocking  force.  In  this  case  the  adaptive  systems  are  near  independent  as  the  maximum 
disturbance  by  neighbouring  actuations  was  less  than  1  degree. 


Table  6  Unit  cells  deflection  under  sequential  actuation 


Cell  1  (deg) 

Cell  2  (deg) 

Cell  3  (deg) 

Total  (deg) 

Activation  Stage  1 

3.7 

0 

0 

3.7 

Activation  Stage  2 

3.8 

3.3 

0 

7.1 

Activation  Stage  3 

3.8 

3.3 

2.3 

9.4 

De-activation  Stage  1 

3.8 

3.2 

0 

7.0 

De-activation  Stage  2 

3.7 

0 

0 

3.7 

De-activation  Stage  3 

0 

0 

0 

0 
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The  servo  motor  eontroller  (Pololu  USB- 16)  offers  8-bits  eontrol  eommand  that  ereates  256 
positioning  options  for  eaeh  servo.  As  eoneluded  previously  in  the  previous  seetion  the 
eurrent  adaptive  system  has  almost  independent  aetuation  amongst  eaeh  eell.  Therefore  the 
total  eombination  of  eamber  profiles  that  ean  be  aehieved  by  the  eurrent  3 -eell  adaptive 
system  is 

256^=16.78^6  (7) 

Of  eourse  this  number  is  heavily  dependent  upon  the  number  of  eells  in  the  adaptive  system, 
and  also  the  resolution  of  eaeh  independent  aetuation  system.  Figure  17  illustrates  the 
adaptation  range  of  this  3 -eell  adaptive  eamber  system,  the  region  enelosed  by  the  upper  and 
lower  eurve  is  the  region  of  operation,  where  as  the  regions  outside  the  eurve  is  the  eurrently 
un-attainable.  As  indieated  by  the  plot  this  partieular  adaptive  system  ean  aehieve  a  eamber 
range  of  +  9.4  degrees  at  the  trailing  edge. 


Figure  17  Adaptation  range  of  the  3-cell  adaptive  camber  system 


3.4  Summary 

From  manufacturing  and  testing  of  the  prototype,  key  concepts  envisioned  by  the  “unit  cell 
approach”  were  validated.  First  of  all,  the  core  concept  of  using  unit  cell  structure  to  achieve 


27 


large  adaptation  through  the  accumulation  of  individual  unit  cell  adaptation  was  validated. 
The  current  3 -cell  adaptive  system  is  capable  of  reaching  total  deflection  of  +9.4  deg  at  the 
trailing  edge  by  accumulating  the  deflection  angle  achieved  by  each  individual  cell.  The 
implication  of  this  achievement  is  far  reaching;  in  an  adaptive  system  where  the  adaptation  is 
distributed  over  a  wide  area  the  shape  change  is  insured  to  be  smoother  and  thus  for  an  aero- 
vehicle  more  aerodynamically  favourable  than  an  adaptive  system  with  lumped  actuation.  A 
further  advantage  of  this  is  that  it  reduces  the  stress  and  strain  concentration  of  the  deformed 
region,  reduces  the  energy  requirement  on  the  actuator  and  widens  the  possible  candidates 
for  both  material  and  actuator  selection.  Secondly,  the  aim  of  achieving  independent,  un¬ 
coupled  individual  adaptive  control  was  realized  in  this  design.  The  actuation  of  each  unit 
cell  has  no  impact  on  the  performance  of  other  unit  cells  and  the  total  actuation  performance 
is  the  linear  combination  of  the  individual  adaptive  units.  This  feature  needs  to  be 
implemented  in  the  actuator-structure  integration  by  ensuring  the  actuators  provide  sufficient 
blocking  force.  In  this  particular  case  the  servo  motors  provided  enough  blocking  force  to 
nullify  the  interaction  between  actuators,  with  maximum  disturbance  in  performance  of  1 
degree.  The  benefit  of  having  multiple  sources  of  independent  controls  is  that  it  enhances  the 
regulation  of  aerodynamic  forces  over  the  aero-vehicle  resulting  in  improved  agility  and 
manoeuvrability.  Thirdly,  the  aim  of  achieving  complex  shape  control  under  high  fidelity 
with  the  unit  cell  design  approach  was  accomplished  with  the  current  design.  The  current 
design  offers  more  than  16  million  unique  shape  profiles  in  just  a  three  cell  network. 
Combined  with  the  previous  point  this  can  be  seen  as  a  possibility  in  moving  away  from 
conventional  fix  wing  design  with  a  few  rigid  control  surfaces  with  lump  control  and  into 
bio-mimic  (bird  flight)  design  with  a  large  number  of  distributed  or  even  continuously  linked 
independent  control  surfaces  that  generate  smooth,  continuous  shape  change  with  large 
numbers  of  degrees  of  freedom. 


4  EVOLUTIONARY  STRUCTURAL  OPTIMIZATION 

4.1  Overview 

It  is  realized  that  the  case  study  conducted  in  section  2  requires  substantial  CPU  time  on  the 
everyday  work  station,  while  at  the  same  time  to  further  perfect  the  design  of  adaptive 
airfoils  additional  schemes,  constraints,  and  optimization  schemes  are  most  likely  required. 
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Hence  at  this  stage  of  the  research  the  computation  cost  of  the  optimization  process  is 
becoming  a  serious  issue,  and  should  be  addressed  appropriately. 

This  section  presents  a  study  that  combines  the  advantages  of  both  the  level  set-based 
boundary  representation  scheme  and  the  ESO  method  to  create  a  new  topology  optimization 
method  named  “Evolutionary  Level  Set’  or  ELS.  The  implicit  dynamic  boundary  scheme  of 
dynamic  level  sets  is  used  to  describe  the  structural  design  boundary,  and  is  coupled  with 
ESO  method  to  advance  the  evolution  of  the  level  set  surface  where  the  design  boundary  is 
embedded  in.  Under  this  scheme,  the  design  optimization  is  changed  in  to  a  numerical 
process  of  updating  the  design  boundary  using  the  evolutionary  method  rather  than  solving 
the  Hamilton-Jacobi  PDEs  via  explicit  time-marching  schemes  (Sethian  1999;  Osher  and 
Fedkiw  2002).  Therefore  the  complex  shape  derivative  analysis  as  well  as  the  numerical 
difficulties  in  most  conventional  level  set  methods  will  be  avoided.  At  the  same  time,  the 
major  advantages  of  ESO  method  such  as  conceptual  simplicity  and  practical  easiness  are 
retained.  It  is  realized  that  the  rationale  for  material  re-distribution  in  ELS  is  largely  based  on 
ESO,  which  as  of  yet  does  not  have  a  comprehensive  theoretical  background.  This  has  been 
subjected  to  critical  reviews  (Querin  and  Rozvany  2002;  Rozvany  2009).  However  in 
practical  engineering  not  being  able  to  fully  understand  a  mechanism  is  no  reason  to  cease 
the  study  and  development  of  such. 

LSM  (Osher  and  Sethian  1988)  was  originally  intended  to  track,  model  and  simulate  the 
evolution  of  dynamic  boundary  with  shape  fidelity  and  topological  change  (Sethian  1 999).  It 
emerged  recently  as  an  alternative  method  to  perform  structural  shape  and  topology 
optimization  without  relaxation  (Sethian  and  Wiegmann  2000).  The  key  concept  underlining 
LSM  is  the  implicit  free  boundary  representation  scheme  by  which  the  design  boundary  is 
described  as  the  zero  level  set  of  a  higher  dimensional  level  set  function.  The  embedded 
boundary  is  mathematically  established  as  a  level  set  model  of  the  Hamilton-Jacobi  partial 
differential  equation  (PDE),  and  a  generic  velocity  field  developed  from  shape  derivative 
analysis  is  included  in  the  PDE  to  enable  dynamic  boundary  advancement.  In  most  LSM 
numerical  difficulties  regarding  the  computational  implementation  of  complicated  Hamilton- 
Jacobi  PDE  need  to  be  taken  into  account  as  summarized  by  (Luo,  Tong  et  al.  2007). 
Furthermore,  a  large  effort  is  devoted  to  deriving  the  design  sensitivity  via  the  shape 
derivative  method.  However,  it  should  be  noted  that  the  implicit  moving  boundary 
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representation  of  level  set  models  itself  is  very  promising  in  representing  complex  interfaces 
by  tracking,  modelling,  and  simulating  the  evolution  of  moving  boundaries  (Luo,  Luo  et  al. 

2006) .  This  representation  has  the  capability  to  simultaneously  achieve  high  shape  fidelity 
and  topological  variation,  and  it  does  not  suffer  from  any  problems  of  explicit  parametric 
surfaces  such  as  wavy  shapes  and  re-parameterization. 

The  ESO  method  (Xie  and  Steven  1993)  differs  from  the  other  methods  in  that  it  is  a 
heuristic  strategy  based  on  the  intuitive  bio-evolutionary  principle,  this  bio-evolutionary 
principle  is  a  natural  phenomenon  of  selecting  the  best  or  fittest  candidate  by  discarding  a 
portion  of  the  least  performing  candidates  of  the  set  in  each  “evolution”  cycle.  ESO  in 
topology  optimization  uses  this  concept  to  structural  analysis  and  achieves  structural 
optimization  by  slowly  removing  inefficient  use  of  material  at  each  iteration.  Most  other 
optimization  methods  employ  a  certain  theoretical  strategy,  such  as  the  OC  (Zhou  and 
Rozvany  1991;  Zhou  and  Rozvany  1992)  and  MMA  (Svanberg  1987;  Svanberg  and  Werme 

2007)  and  search  for  the  optimal  design  point  based  on  gradient  based  sensitivity  analysis 
which  in  theory  insures  the  discovery  of  the  global  optimal.  However,  past  presentations 
such  as  (Querin,  Steven  et  al.  1998)  and  (Huang  and  Xie  2007)  have  shown  that  even  though 
ESO  method  does  not  have  solid  theoretical  bases,  when  applied  appropriately  its  structural 
results  can  be  very  promising  if  not  almost  identical  to  its  counterpart’s  especially  in  load 
bearing  designs.  In  addition,  its  conceptual  simplicity  and  computational  effectiveness 
become  distinctly  attractive  in  scenarios  where  explicit  design  sensitivity  analysis  is  difficult 
or  time  costly  (Nguyen,  Tong  et  al.  2007).  Hence  proving  that  the  field  of  ESO  is  worthy  of 
further  research  and  investigations.  BESO  or  Bi-directional  ESO  (Querin,  Steven  et  al.  1998; 
Li,  Steven  et  al.  2000)  is  an  extension  of  ESO  that  enables  material  addition  as  well  as 
removal,  it  augments  certain  short  comings  in  the  original  “hard  kill”  ESO  such  as  its  initial 
material  distribution  condition,  evolution  stability,  and  optimality. 

4.2  Evolutionary  optimization  philosophy 

Before  EES  is  introduced  it  is  necessary  to  understand  the  core  concept  of  ESO  that  was 
originally  developed.  The  fundamental  concept  of  ESO  (Xie  and  Steven  1993)  is  to  produce 
a  fully  stressed  type  structure  that  has  maximum  stiffness  and  minimum  weight  by  gradually 
removing  local  areas  (finite  elements)  that  are  of  low  utilization  from  the  domain,  i.e. 
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elements  that  are  under  relatively  lower  stress.  The  removal  ean  be  based  on  a  number  of 
numerieal  eonditions,  often  ealled  “rejeetion  eriteria”,  whieh  is  a  performanee  assessment  of 
the  individual  elements  based  on  eertain  iso-value  parameter  that  gives  some  indieation  of 
struetural  performanee.  In  the  past  the  Von  Mises  stress  gvm  and  the  strain  energy  density  H 
have  been  widely  used  (Xie  and  Steven  1997).  Elements  are  to  be  removed  if  they  satisfy  the 
following  eondition 

s[(yvM,e)^S{(y„)  e  =  l:7V  «  =  l:«max 

or  (8) 

g{H^)<g{H„) 

where  ‘‘g”  is  the  funetion  of  the  rejeetion  eriterion,  N  is  the  number  of  elements  in  the 
domain,  nmax  is  the  maximum  iteration  or  loop  eount  of  the  optimization  proeedure,  On  and 
Hnare  the  threshold  values  for  Von  Mises  stress  and  strain  energy  density  at  the  n^^  iteration. 
The  removed  elements  are  to  be  treated  as  voids  and  eease  to  have  further  eontribution  to  the 
struetural  stiffness  by  setting  the  thiekness  te  of  those  elements  to  a  very  small  number.  This 
algorithm  is  repeated  until  eertain  eonstraints  are  met,  sueh  as  equality  eonstraints  and 
eonvergenee. 

Originally  ESO  had  a  “hard  kill”  nature,  in  that  the  finite  element  properties  are  diserete 
(either  solid  or  void)  and  are  irreversible.  This  plaeed  eertain  eonditions  on  ESO  in  order  to 
operate  eorreetly;  sueh  as  in  the  initialization  of  the  EE  domain  areas  that  are  antieipated  to 
be  part  of  the  final  strueture  must  not  be  assigned  as  void  areas;  whieh  ean  be  diffieult  to 
determine.  Or  the  amount  of  material  removed  at  eaeh  iteration  needs  to  be  small  to  avoid 
“over-removal”,  eausing  the  optimal  layout  to  be  missed.  It  is  obvious  what  is  desired  is  to 
let  the  finite  element  properties  to  be  reversible,  (Querin,  Steven  et  al.  1998;  Yang,  Xie  et  al. 
1999)  proposed  a  BESO  method  whieh  adds  additional  eriteria  to  the  performanee 
assessment  to  allow  void  areas  to  turn  into  solids,  making  the  optimization  type  still  diserete, 
but  reversible.  This  method  works  well  when  eertain  pre-requisite  eonditions  on  the 
initialization  of  the  domain  are  satisfied,  and  eurrent  development  is  underway  to  resolve  the 
sensitivity  issue  surrounding  the  employed  rejeetion  eriteria.  (Li,  Steven  et  al.  2000)  made 
modifieations  to  the  original  ESO  to  make  the  optimization  partially  eontinuous  by  step¬ 
sizing  te  between  the  bound  of  tmin  and  tmax,  and  used  the  sensitivity  analysis  of  the  objeetive 
funetion  with  respeet  to  the  thiekness  as  part  of  the  rejeetion  eriteria.  The  modified  algorithm 
expanded  the  ESO  realm  in  to  more  eomplex  design  seenarios  by  making  appropriate 
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tradeoffs  in  computation  complexity  and  CPU  cost.  Other  forms  of  bidirectional  ESO  such 
as  (Querin,  Young  et  al.  2000;  Liu,  Jin  et  al.  2008)  are  generally  based  on  the  above  two. 

4.3  Implicit  boundary  representation  via  level  sets 

As  aforementioned,  the  key  concept  of  the  level  set  models  is  to  develop  an  implicit  free 
boundary  which  is  mathematically  described  by  a  first-order  Hamilton- Jacobi  equation 
(Sethian  1999) 

^^^  +  F|V(1)|  =  0,  cl)(x,0)  =  cl)o(x)  (9) 

ot 

where  0(x,^)  is  the  Lipschitz  continuous  scalar  function  and  t  is  the  pseudo-time  to  enable 
the  dynamic  process  of  the  level  set  surface,  and  V  is  the  normal  velocity  field.  In  this  work, 
the  movement  of  the  boundary  is  achieved  by  transporting  the  level  set  function  according  to 
structural  strain  energy  density  rather  than  a  series  of  solutions  from  the  Hamilton- Jacobi 
PDE  (Wang,  Wang  et  al.  2003).  Hence,  this  study  does  not  need  to  solve  the  Hamilton- 
Jacobi  PDE  using  explicit  schemes  and  calculate  the  velocity  V  using  shape  derivative 
method  (Allaire,  Jouve  et  al.  2004).  Instead,  the  evolutionary  algorithm  based  on  the  concept 
of  ESO  is  applied  to  update  the  level  set  surface  O  as  well  as  the  design  boundary  in  a 
similar  manner  as  a  fully  stressed  structural  design. 

Hence,  a  major  attribute  of  EES  is  the  combination  of  the  level-set  boundary  representation 
scheme  and  the  bio-evolutionary  principle,  leading  to  a  smooth,  distinct  structural  design 
boundary.  In  this  method  all  nodal  points  of  the  finite  element  domain  are  evaluated  by  the 
function  O  that  gives  an  indication  of  the  relative  structural  performance.  A  higher  value 
indicates  that  this  part  of  the  structure  has  high  performance,  while  a  low  value  would 
indicate  that  the  region  is  of  low  performance.  A  3D  level  set  contour  can  be  constructed 
with  the  O  function  over  the  entire  domain  to  give  indication  of  relative  performance  of  the 
local  design  domain  (see  Figure  18).  An  iso-surface  Sn  can  then  be  determined  to  intercept 
the  level  set  contour  at  a  chosen  magnitude  at  the  n^^  iteration,  thus  causing  level  set  contours 
of  the  high  performance  regions  to  lie  above  the  iso-surface  while  the  level  set  contours  of 
the  low  performance  regions  to  be  below.  The  structural  topology  can  then  be  distinctly 
determined  by  the  2D  mapping  of  the  level  set  contours  above  the  iso-surface.  This  ESM 
inspired  implicit  boundary  representation  scheme  generates  a  smooth,  distinct  boundary  that 
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clearly  separates  solid  and  void  regions.  Whereas  most  current  topology  optimization 
methods  can  achieve  distinct  boundaries,  but  lack  the  smooth  surface  or  require  post¬ 
processing. 


ELS’  material  distribution  scheme  follows  the  concept  of  classic  ESO  (Xie  and  Steven  1997) 
in  removing  material  at  lowly  utilized  areas,  but  as  shall  be  explained  in  detail  in  following 
sections  ELS’  unique  bi-directional  algorithm  also  distributes  material  based  on  its 
performance.  The  material  distribution  of  ELS  is  based  on  element  level  but  assessed  by  the 
nodal  performance  of  the  element.  Individual  finite  elements  with  all  nodal  O  functions 
greater  than  the  iso-surface  Sn  are  assigned  a  stiffness  weighting  factor  of  unity  to  represent 
solid,  finite  elements  with  all  nodal  <I)  functions  below  Sn  are  assigned  a  factor  of  zero. 
However,  elements  with  nodal  O  functions  both  above  and  below  the  iso-surface  Sn  (i.e. 
boundary  elements)  are  assigned  a  weighting  based  on  the  “Ersatz”  model.  The  use  of 
“Ersatz”  model  (Allaire,  Jouve  et  al.  2004)  effectively  makes  ELS  a  continuous  optimization 
method,  the  first  from  the  ESO  branch  according  to  available  publications.  Being  a 
continuous  optimization  method  carries  many  benefits  besides  allowing  implicit 
representation  of  boundary,  it  carries  the  benefit  of  avoiding  design  variable  explosion  that 
some  discrete  methods  have,  plus  it  allows  the  possibility  to  apply  theoretically  well-founded 
gradient-based  optimization  algorithms  such  as  the  OC  (Zhou  and  Rozvany  1991)  and  MM  A 
(Svanberg  1987). 
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4.4  Problem  formulation  with  ELS 

The  optimization  of  load  bearing  static  structures  by  most  topology  optimization  methods  are 
referred  to  as  mean  compliance  problems,  because  the  objective  function  -  the  mean 
compliance  is  being  minimized.  By  principle  of  virtual  work  under  equilibrium  condition  the 
mean  compliance  is  equivalent  to  the  internal  virtual  work.  Therefore  for  an  elastic  body 
with  arbitrary  shape  this  is  expressed  as 

+  ^tiSudT 

n  nr 

Where  f  and  f  are  the  vectors  of  body  force  and  tractions,  imposed  on  region  Q  and  side  T, 
while  u  is  the  actual  displacement  field,  8u  is  the  virtual  displacement  field.  Dpi  is  the 
stiffness  tensor,  ay  and  Ski  are  the  strain  and  virtual  strains  respectively. 

As  for  ESO  and  ELS  method,  the  optimal  state  is  different  to  the  above  and  is  defined  by  a 
fully  stressed  type  state  rather  than  a  minimal  compliance  state,  and  for  some  ESO  methods, 
a  so  called  “performance  index”(Querin,  Steven  et  al.  1998).  The  validity  of  these  methods 
have  been  examined  in  the  past  (Rozvany  and  Querin  2004;  Chiandussi  2006;  Rozvany 
2009)  ,  sometimes  critically.  However  numerous  literatures  such  as  (Tanskanen  2002)  and 
(Patnaik  and  Hopkins  1998)  have  also  argued  with  logical  theoretical  bases  for  the 
optimality  of  the  ESO  method  and  the  fully  stressed  design  (FSD)  approach.  Regardless,  it  is 
not  the  aim  of  this  report  to  take  sides  in  the  optimality  debate  but  to  present  a  new  topology 
optimization  method,  which  is  given  below. 


The  generic  problem  statement  for  ELS  is  shown  below 

Max:  min(0(x))  Vx=(x,. ,  y, )  e  50*  n  Q* 

Tf  =  KU  Or  equilibrium  equantions 

r  -  (10) 

Subject  to:)  J^.  w^dLl  =  V  i  =  1,2, 3. ...n 

or  0(x)  >  O 

The  objective  function  is  max(min(0(x))) .  The  function  <1)  evaluates  to  an  iso-value  that  can 
be  used  to  compare  the  performance  of  each  node,  potential  O  functions  in  load  bearing 
designs  are  strain  energy  density,  and  Von  Mises  stress.  0(x)  is  the  d)  function  value  of  all 

points  in  x,  for  x=(X;,y,  )  g  50*  nO*  which  represents  the  solid  regions  of  the  design  domain. 
To  maximize  the  minimum  d)  function  in  the  solid  domain  is  to  lift  the  performance  of  the 
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structural  region,  the  optimal  solution  is  when  each  component  of  the  structure  has  achieved 
its  full  performance  potential.  The  topology  of  the  structure  is  expressed  implicitly  as 
follows  and  is  illustrated  in  Figure  19 

solid  Q*  =  :  0(x)  > 

boundary  dQ*  =  =  (H) 

void  Q/Q*  =  :  0(x)  < 

In  ELS  the  3D  level  set  contour  is  formed  by  the  O  value  of  each  node  in  the  design  domain, 
consequently  the  active  design  space  is  implicitly  constructed  by  the  2D  mapping  of  the 
partial  level  set  contour  that  satisfies  the  condition  o(x)>^„  on  to  the  iso-surface  Sn,  where 

Sn  is  the  iso-surface  at  the  n-th  iteration  of  the  optimization  process,  {f}  =  [^]{l/}  is  the 

governing  equation  of  the  system,  Wi  is  the  element  weighting  factor,  V  is  the  volume  of 
allowable  material  and  O  is  the  minimal  allowable  parameter. 


Figure  19  Design  domain  of  ELS 

The  subsequent  numerical  procedure  involves  the  finite  analysis  of  the  design  domain.  For  a 
linear  mechanical  system  this  would  be 

=  (12) 

where  [^]  is  the  global  stiffness  matrix,  and  {t/}  and  {F}are  the  global  displacement  and 

force  fields.  When  constructing  the  element  stiffness  matrices  the  element  weighting  factor 
We  of  each  element  is  applied  to  the  integration  as  shown  in  Eq.  (13)  to  distinguish  between 
solid  and  void  region.  The  weighting  factor  acts  similarly  to  the  relative  material  density  p  in 
SIMP  (Bendsoe  and  Sigmund  2003). 

W  ■  <  W  <  1 


35 


As  indicated  by  the  equation,  the  ELS  is  effeetively  a  continuous  optimization  method.  This 
presentation  is  limited  to  linear  elasticity  only  for  the  purpose  of  numerical  simplicity, 
though  there  is  no  conceptual  difficulty  in  the  employment  of  nonlinear  model. 


Since  ELS  evolves  around  the  fundamental  concept  of  maximum  utilization  of  the  load 
bearing  capabilities  of  the  structure,  this  involves  solving  the  <1)  function  that  indicates  the 
effeetive  utilization  of  eaeh  finite  element,  or  more  specifically  at  each  nodal  point  of  the 
finite  element  if  constant  strain  elements  are  not  used;  typically  the  von  Mises  stress  and 
the  strain  energy  density  H  can  be  used  to  compare  the  effective  utilization  of  each  element, 
but  that  does  not  mean  alternative  parameters  can  not  be  used.  The  finite  element  form  of  the 
von  Mises  stress  and  the  strain  energy  density  of  a  quad-4  element  is  given  by  Eq.  (14). 


(14) 


where  e  denotes  the  current  finite  element,  (o-J  andj^J  are  the  stress  matrix  and  strain  matrix 
respectively,  [c]  is  the  von  Mises  coefficient  matrix  given  by  Eq.  (15)  for  an  quadratic 
element. 


[c-] 


1 

-.5 

0 


-.5  0 

1  0 

0  3 


(15) 


The  value  at  the  nodal  points  is  evaluated  by  averaging  the  Gaussian  quadrature  in  the 
numerical  integration  using  the  Gaussian  points  surrounding  the  node  as  shown  in  Figure  20. 

Evaluated  using  surrounding  Gaussian  points 


Gaussian  point 
Nodal  point 


Figure  20  Evaluation  of  nodal  attributes  through  surrounding  Gaussian  points 
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For  example  to  interpolate  the  nodal  value  rather  than  the  element  value,  the  nodal  stress  at 
node  “j”  |(TyJ  is  evaluated  by  the  averaged  Gaussian  quadrature.  Therefore  instead  of  basing 

the  construction  of  the  strain  displacement  matrix  [B]  on  the  four  Gaussian  points  within  the 
element,  the  new  [B]  is  based  on  the  Gaussian  points  around  the  node  as  in; 


E,(a)’^EAa) 


Ej(a)  =  {ei,e2,...eJ 


(16) 


where  m  is  the  number  of  Gaussian  points  around  the  node,  Ej  contains  the  set  of  elements 
that  share  the  node  j,  ^  and  i]  are  the  natural  coordinates  of  the  Gaussian  point  from  each 
element  in  E/a). 


By  averaging  the  Gaussian  quadrature,  combined  with  the  continuous  effect  of  the  element 
weighting  factor,  the  optimization  method  has  the  ability  to  revert  void  regions  back  to  solid 
or  partially  solid  which  eliminates  the  adverse  effect  from  over-removal  and  reduces  the 
possibility  of  sub-optimal  solutions.  More  on  this  issue  is  explained  in  section  5. 

4.5  Numerical  implementations 

Creating  iso-surface 

From  the  graphical  illustration  in  Figure  18  the  iso-surface  Sn  is  a  transverse  plane  in  x  and  y 
raised  to  a  certain  height  in  z.  It  acts  as  the  one  and  only  rejection  criterion  required  for  ELS. 
Finite  elements  associated  with  O  values  which  lie  completely  below  the  iso-surface  are  to 
have  their  element  weighting  factor  zeroed,  which  effectively  turns  the  region  into  void. 
Finite  elements  associated  with  O  values  that  lie  partially  below  the  iso-surface  (at  least  one 
nodal  attribute  above  the  iso-surface)  are  to  have  the  element  weighting  factor  determined 
using  the  Ersatz  model  explained  later  in  the  chapter.  Elements  associated  with  O  values  that 
lie  above  Sn  are  treated  as  solid  and  have  weighting  factor  valued  as  1 . 


For  problems  with  volume  equality  constraint,  the  rejection  criterion  can  be  imposed  in  a 
fixed  fashion  or  a  relaxed  fashion.  The  fixed  rejection  criterion  is  established  by  setting  iso¬ 
surface  Sn  to  implicitly  impose  volume  equality  constraint,  this  can  be  done  by  the  following 
procedure. 

Step  1.  Create  a  matrix  of  N  by  2,  with  one  column  being  the  nodal  iso-value  Oj  and  the 
other  the  nodal  density  wj 
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Step  2.  Sort  the  matrix  according  to  the  Oj  column  in  descending  fashion;  this  can  he  done  in 
MATLAB  effectively  using  the  SORTROW  function 

Step  3.  Sum  the  nodal  density  wj  starting  from  the  first  row,  and  Sn  is  the  Oj  value 
corresponding  to  the  row  at  which  the  summation  of  the  nodal  density  wj  first  violates  the 
volume  equality 

The  relaxed  rejection  criterion  is  simply  hasing  the  selection  of  Sn  on  the  volume  difference 
between  the  starting  material  volume  and  the  constraint,  so  that  at  each  iteration  it  reduces 
the  volume  difference  with  the  eventual  goal  of  reducing  the  volume  difference  to  zero  in  a 
certain  number  of  iterations.  The  reduction  of  the  volume  difference  can  be  done  in  fixed 
increments  or  in  a  continuous  manner.  Regardless  of  using  fixed  or  relaxed  rejection 
criterion,  the  final  topologies  are  identical.  However,  the  relaxed  criterion  converges  faster  in 
some  cases.  For  problems  with  stress  or  strain  energy  density  equality  constraint,  the  fixed 
rejection  criterion  is  established  by  setting  Sn  to  the  equality  constraint  at  the  value  O  ,  while 
the  relaxed  criterion  can  be  done  similar  to  that  of  the  volume  condition. 


Reassessing  weighting  factor 

The  expression  for  the  elastic  field  of  the  current  state  with  respect  to  the  intrinsic  material 
property  of  the  structure  is  represented  using  the  well  recognized  Ersatz  material  model.  In 
the  Ersatz  material  model,  elements  with  all  nodes  located  below  the  iso-surface  are  all 
deemed  to  be  void  occupying  regions  implemented  by  mimicking  artificially  weak  phases 
without  risking  numerical  singularity,  elements  with  all  nodes  located  above  the  current  iso¬ 
surface  are  treated  as  solid  with  the  respective  relative  density  as  unity.  The  boundary 
element’s  (with  a  mixture  of  nodes  above  and  below  the  iso-surface)  weighting  factors  are 
interpolated  based  on  the  area  ratio  of  the  2D  projection  of  the  plane  which  formed  by  the  O 
values  of  the  nodal  points  that  lies  above  the  iso-surface  and  the  total  area  of  the  element,  as 
shown  in  Figure  21.  The  element  weighting  factors  for  such  boundary  elements  can  be  easily 
calculated  as 


A  =  A  -A~ 


(17) 


Where  A^  and  A  are  the  projected  2-D  areas  of  the  level  set  surface  that  lies  above  and 
below  the  iso-surface,  and  Ae  is  the  total  area  of  the  finite  element. 
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I 


Iso-surface 


Figure  21  Ersatz  model  for  interpolating  weighting  factor  of  boundary  elements 


After  the  establishment  of  a  new  iso-surface  and  the  interpolation  of  the  relative  densities,  a 
dual  check  is  made  to  see  first  whether  the  constraint  is  met  and  second  if  the  solution  has 
converged.  If  not  the  procedure  goes  through  another  round  of  finite  element  analysis.  The 
convergence  criterion  is  non-unique,  the  one  used  in  this  report  is  given  in  equation  (18). 


\s-s^ 


-  <  error 


(18) 


Filtering 

Likewise  to  other  topology  optimization  procedures  based  on  finite  element  methods  the 
ELS  is  prone  to  checkerboard  features  (which  is  regarded  as  non-manufacturable  solution) 
but  to  a  significantly  less  degree.  Also  likewise  to  other  topology  optimization  procedures, 
this  issue  can  be  suppressed  by  simple  filtering  techniques.  The  filtering  techniques  are 
numerous  but  fundamentally  they  achieve  the  same  purpose  of  blurring  or  smoothing  the 
weighted  matrix  or  the  sensitivity  of  the  weighted  matrix.  In  this  work  a  convolution  kernel 
is  used  to  filter  the  element  weighting  factors.  Its  effectiveness  is  illustrated  below  in  Figure 
22.  It  is  interesting  to  note  that  even  in  the  un-filtered  topology  the  checkerboard  features  are 
at  a  minimal  compared  to  other  topology  optimization  methods,  this  is  because  in  the  process 
of  constructing  the  level  set  surface  <!>  the  nodal  parameter  was  evaluated  by  the  Gaussian 
points  surrounding  the  node  i.e.  from  the  neighbouring  elements,  thus  in  a  way  it  has  already 
achieved  a  “smoothing”  effect. 
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Figure  22  Result  of  un-filtered  and  filtered  topologies  for  tip  loaded  cantilever 

Topology  growth  and  removal 

Fundamentally  the  removal  philosophy  in  ELS  is  similar  to  the  approach  of  the  classic  ESO 
method  which  is  to  remove  the  non-performing  or  low  performing  regions.  However,  in  this 
method  local  material  growth  and  not  just  material  removal  also  occurs  during  the 
optimization,  which  is  due  to  the  combined  effect  of  the  continuous  element  weighting  factor 
and  the  averaged  Gaussian  evaluation.  For  example  in  Figure  23  where  the  weighting  factors 
of  element  A,  B  and  C  are  all  non-zero,  but  zero  for  D.  In  situations  where  the  nodal  attribute 
of  any  of  the  three  nodes  marked  by  the  dotted  circle  are  above  the  current  iso-surface  Si, 
then  by  Ersatz  method  the  weighting  factor  of  element  D  will  increase  to  a  non-zero  value  at 
the  next  iteration,  allowing  the  structure  to  ‘‘grow”,  thus  creating  the  effect  of  removing  the 
under  performing  regions  and  strengthens  the  high  performing  regions.  Figure  24  is  an 
example  of  a  cantilever;  the  figures  show  the  structure  evolving  by  undertaking  local 
material  growth  and  removal  while  maintaining  a  constant  overall  volume. 
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Figure  23  Illustration  of  element  growth 
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step  i-1  Step  / 

Figure  24  Material  re-distribution  process  mapped  using  W 


4.5  Examples 

In  this  section  3  numerical  examples  are  given,  and  each  is  presented  with  complete  list  of 
results  and  discussions.  The  examples  are  classic  examples  of  MBB  and  cantilever  beam,  for 
which  the  optimized  topologies  are  well  known  and  can  be  compared  with  the  ELS  result  to 
validate  its  accuracy  and  practicality. 

As  mentioned  in  the  previous  sections,  the  ELS  method  offers  greater  robustness  in  the 
initialization  of  design  domain.  This  is  showcased  in  the  examples  given;  all  examples  are 
initialized  with  random  values  of  material  distribution  unless  otherwise  stated.  The  domains 
are  meshed  using  1mm  x  1mm  quad4  elements  in  all  3  cases.  The  O  function  in  this  chapter 
is  the  strain  energy  density.  The  material  used  has  an  elastic  modulus  of  lOOOMPa,  with 
Poisson  ration  of  0.3 1 
Case  1 

This  example  is  a  cantilever  with  tip  load  of  100  N  placed  on  its  neutral  axis.  The  domain  is 
of  size  120mm  x  40mm.  The  imposed  constraint  on  the  optimization  is  a  volume  equality 
constraint  of  50%,  this  is  imposed  as  fixed  rather  than  relaxed.  The  problem  is  illustrated  in 
Figure  25.  The  final  topology  as  well  as  some  of  the  intermediate  solutions  are  shown  in 
Figure  26,  their  corresponding  level  set  contours  are  shown  in  Figure  27.  The  plot  of  iso¬ 
surface  history  is  shown  in  Figure  28.  The  total  strain  energy  history  is  shown  in  Figure  29. 
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Figure  26  Topology  history  for  case  1  and  final  iteration 
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Figure  27  Level  set  contour  of  strain  energy  density  distribution 


Figure  28  Plot  of  iso-surface  history 


Figure  29  Plot  of  total  strain  energy  history 


The  topology  starts  out  as  a  simple  beam  structure  with  a  single  truss  network  near  its 
boundary  supports.  As  the  optimization  progresses  additional  truss  networks  gradually 
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appear,  while  the  top  and  bottom  layer  of  the  beam  redistribute  the  material  away  from  the 
tip  and  towards  the  wall.  The  final  total  strain  energy  of  the  structure  resolved  to  be  1.82  x 
10^  mJ.  The  strain  energy  is  distributed  evenly  throughout  most  of  the  structure,  as  indicated 
in  Figure  27,  but  increases  shapely  near  the  boundary  supports  and  the  loading  point,  the 
increase  in  strain  energy  can  be  order  of  one  magnitude. 

Case  2 

The  second  case  is  a  cantilever  with  a  tip  load  of  50  N.  The  imposed  constraint  on  the 
optimization  is  a  strain  energy  density  inequality  constraint  H  >  0.3  mJ/mm^.  This  is 
imposed  as  a  fixed  constraint.  The  problem  is  illustrated  in  Figure  30.  The  resulting  topology 
should  tell  the  designer  the  amount  of  volume  that  is  required  to  meet  the  strain  energy 
density  constraint. 


Figure  31  displays  the  final  topology  and  some  intermediate  solutions,  and  their 
corresponding  strain  energy  density  contour  is  shown  in  Figure  32.  Figure  33  shows  the 
volume  fraction  history.  Figure  34  displays  the  total  strain  energy  history. 
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Figure  31  Topology  history  for  case  2  and  final  iteration 


Figure  32  Level  set  contour  of  final  strain  energy  distribution 


Figure  33  Plot  of  volume  fraction  history 
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Figure  34  Plot  of  total  strain  energy  history 


The  figures  show  that  ELS  converged  to  the  well  expected  topology.  The  volume  fraction 
ended  up  being  at  40%  after  more  than  200  iterations.  This  is  due  to  the  convergence 
criterion  not  being  met  (perhaps  too  demanding),  even  though  the  topology  changes  are 
hardly  detectable  after  120  iterations  or  so.  The  final  total  strain  energy  of  the  structure 
resolved  to  be  0.58  x  10^  mJ. 

Case  3 

The  third  example  is  a  classic  case  of  MBB  design,  with  a  centre  downward  force  of  50  N. 
The  imposed  constraint  on  the  optimization  is  a  volume  equality  constraint  of  25%  of  total 
design  domain.  In  this  case  the  equality  constraint  is  relaxed  and  the  domain  is  initialized  as 
full.  The  design  domain  is  illustrated  in  Figure  35. 


Figure  35  Design  domain  for  case  3 

Figure  36  and  Figure  37  displays  the  topology  history  and  the  final  strain  energy  contour. 
Because  the  volume  equality  constraint  is  imposed  in  a  relaxed  fashion,  the  volume  ratio  of 
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the  design  gradually  changes.  The  history  of  volume  fraction  change  is  shown  in  Figure  38, 
and  the  iso-surface  history  is  shown  in  Figure  39,  its  total  strain  energy  plot  is  in  Figure  40 


Figure  36  Topology  history  for  1  (top  left),  20  (top  right),  50  (bottom  left)  and  final  iteration 


Figure  37  Level  set  contour  of  final  strain  energy  distribution 


Figure  38  Volume  fraction  history 
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iteration 

Figure  40  Plot  of  total  strain  energy  history 

The  topology  history  shows  an  interesting  pattern,  with  the  initial  stages  representing  a  main 
support  arc  with  numerous  bars  transferring  the  load  to  the  arc.  However,  as  the  constraint  is 
continuously  “tightened”,  the  arc  straightens  and  eventually  becomes  a  trapezium  shaped 
truss  network.  The  final  total  strain  energy  of  the  structure  resolved  to  be  0.23  x  10^  mJ. 

The  results  presented  in  all  three  cases  concurs  with  other  pre-existing  optimization  methods 
both  topologically  and  numerically,  minor  differences  are  most  likely  due  to  different 
filtering  techniques  applied.  It  is  impractical  to  accurately  compare  the  relative 
computational  efficiency  of  different  optimization  methods  in  all  aspects,  but  to  compare  the 
CPU  run  time  in  general  terms  a  modest  description  of  ELS  method  would  be  “competitive” 
to  the  other  optimization  methods  using  the  same  finite  element  mesh. 

4. 6  Summary 

This  section  has  presented  a  new  topology  optimization  method  named  “ELS”,  which 
follows  the  major  concept  of  ESO  philosophy  but  with  its  unique  formulation  and  described 
with  an  implicit  free  boundary  representation  scheme  of  dynamic  level  sets.  The  unique 
algorithm  allows  material  distribution  (both  removal  and  reinforcement)  to  be  determined 
based  on  a  simple  but  efficient  criterion,  with  the  implicit  effect  of  removing  material  from 
low  efficiency  regions  and  adding  material  around  high  efficiency  regions.  The  resulting 
topologies  from  ELS  are  distinctly  marked  in  solid  and  void  and  have  smooth  boundaries. 
They  also  indicate  that  the  final  structure  is  in  a  state  of  nearing  “fully  stressed”  saved  for  the 
boundary  support  regions  and  the  loading  point.  The  examples  show  ELS  to  be  effective  in 
generating  optimal  topologies  at  a  reasonable  convergence  rate. 
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5.  CONCLUDING  REMARKS 


The  main  features  of  this  report  include  multi-objective  optimization  of  tapered  unit  cell 
airfoil  structures  for  camber  deflection,  construction  and  testing  of  a  prototype  unit  cell 
network,  and  the  introduction  of  a  new  topology  optimization  algorithm  referred  to  as 
“Evolutionary  level  set”  for  the  role  of  optimizing  load  bearing  continuum  structures. 

The  method  of  conducting  multi-objective  optimization  on  the  tapered  unit  cell  structure  was 
proven  by  comparison  to  be  superior  than  single  objective  optimization  in  delivering 
practical  designs  that  would  generate  sufficient  adaptation,  while  exhibit  adequate  stiffness, 
ft  is  also  found  that  in  general  stiffness  favoured  designs  are  more  “design  friendly” 
compared  to  adaptation  favoured  designs  by  manifesting  relatively  less  complex  truss 
structures  with  the  majority  of  structural  members  taking  axial  loads,  and  in  addition  these 
designs  tend  to  feature  relatively  simple  actuation  mechanisms  which  are  easy  to  implement 
in  the  physical  world.  Where  as  for  the  adaptation  favoured  designs  the  main  problem  with 
actuation  was  in  fulfilling  the  uni-morph  mechanism  in  practice. 

The  feasibility  of  the  tapered  unit  cell  design  was  demonstrated  and  proven  by  prototype 
testing.  Two  major  findings  arise  from  the  testing.  First  of  all,  the  “generating  sizable 
adaptation  through  accumulation  of  individual  adaptation”  condition  for  the  unit  cell  was 
confirmed  in  the  testing.  The  individual  adaptation  of  each  of  the  three  cells  in  the  network 
was  no  more  than  3.8  degrees,  but  when  actuated  in  unity  the  a  sizeable  deflection  of  9.4 
degrees  that  was  build  up  iteratively  from  each  cell  was  reached  at  the  trailing  edge.  The 
second  finding  is  the  confirmation  of  independent  control  network.  In  the  testing  the  control 
input  issued  to  any  cell  does  not  impede  or  compromise  the  control  capability  of  other  unit 
cells  in  the  network  (Given  that  the  actuators  does  not  suffer  from  un-wanted  deformation 
and  have  sufficient  blocking  force).  This  fact  leads  to  the  total  number  of  configurations  of 
the  adaptive  airfoil  being 

(19) 

Where  At  is  the  number  of  possible  total  configuration,  Ac  is  the  number  configuration  for 
each  actuator,  and  N  is  the  total  number  of  actuators  in  the  system.  It  is  straightforward  to 
see  that  unit  cell  design  offers  a  highly  adaptable  system  with  the  possible  number  of 
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adaptation  configurations  being  exponentially  proportional  with  respect  to  the  number  of 
unit  cells  within  the  network. 

The  new  optimization  algorithm  ELS  exhibits  potential  in  the  design  of  aero-vehicle 
structures  by  presenting  highly  detailed  structural  geometry  at  a  low  computational  cost.  The 
implicit  boundary  representation  scheme  in  Evolutionary  Level  Set  and  other  Level  Set 
methods  does  not  suffer  from  negative  impact  of  blurring  boundaries  (from  using  filtering 
techniques)  and  zigzag  boundaries  (from  using  discrete  “hard  kill”  methods).  Simultaneously 
Evolutionary  Level  Set  drastically  reduces  the  computation  cost  by  avoiding  the  numerically 
intensive  Hamilton-Jacobi  equation  and  implementing  a  less  numerically  intensive 
Evolutionary  Structural  Optimization  algorithm.  The  impact  on  the  result  of  the  optimization 
is  that  the  optimization  is  no-longer  a  mean  compliance  optimization  that  has  strong 
theoretical  proof  of  optimality  and  favoured  by  academic  researchers,  but  rather  the 
optimization  is  a  fully  stress  optimization  that  does  not  yet  has  the  same  level  of  theoretical 
background  as  the  mean  compliance  but  is  favoured  and  deep-rooted  in  the  engineering 
industry. 

In  this  final  section  of  the  report  a  design  framework  is  presented  to  summarise  the  design 
environment  of  the  adaptive  airfoil.  The  proposed  design  framework  is  garnered  from  the 
research,  design,  manufacturing,  and  testing  experience  gained  throughout  this  project.  The 
design  framework  is  formed  by  three  interdependent  design  considerations,  forming  a  triad 
consideration  factor  consisting  practical  consideration,  structural  consideration,  and  material 
consideration  (Ligure  41). 
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Practical  considerations  include  design  matters  that  influence  the  practicality  of  the  design  in 
a  real  world,  commercial  sense.  It  includes  matters  such  as 

•  Cost  requirement  on  the  adaptive  airfoil  in  order  to  be  not  only  financially  viable,  but 
commercially  competitive  against  other  performanee  enhancing  methods 

•  Regulatory  requirement  set  by  the  industry,  government  body,  or  international 
institution  on  restrictions  of  aircraft  related  matters,  such  as  wing  span,  noise  control, 
green  house  gas  emission,  etc 

Structural  considerations  are  assoeiated  with  the  performance  requirement  on  the  wing 
strueture,  ineluding  both  the  passive  strueture  and  the  integrated  actuator.  Examples  are 

•  Adaptation  requirement  in  the  magnitude  of  shape  change  the  structural  system 
should  provide 

•  Stiffness  requirement  in  the  magnitude  of  un-wanted  deformation  the  structural 
system  is  allowed  to  sustain  under  external  loading 

Material  considerations  are  associated  with  the  material  property  requirement  on  both  the 
passive  structure  and  the  actuator;  this  is  different  to  the  structural  considerations  in  that 
struetural  performance  sueh  as  adaptation  and  stiffness  can  be  altered  by  design,  whereas 
material  performance  can  not  be  altered  and  is  an  intrinsic  property  of  the  material. 
Examples  are 
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•  Flexibility  requirement  on  the  structural  materials  in  the  adaptive  system,  this 
includes  both  the  yield  strain  parameter  on  the  host  material  and  the  actuation  limit  of 
the  actuator 

•  Strength  requirement  on  the  structural  materials  in  the  adaptive  system,  which 
includes  the  tensile  stress  value  of  the  passive  material  and  the  blocking  force  of  the 
actuator 

As  aforementioned  the  triad  factors  are  interdependently  coupled,  therefore  there  are  design 
issues  associated  with  not  to  a  single,  but  two  or  all  three  factors.  Examples  of  this  are 

•  Density/weight  requirement  on  the  adaptive  system  is  a  material  consideration  as  the 
material  density  is  an  intrinsic  material  factor,  but  the  weight  of  the  airfoils  account 
for  a  significant  portion  of  the  total  aircraft  weight  (10%  for  passenger  jet)  and  has  a 
sizable  impact  on  the  financial  margins  in  its  commercial  operations.  Thus  it  is  also 
practical  consideration 

•  Energy  efficiency/consumption  on  the  adaptive  system  is  a  requirement  that  overlap 
all  three  factors,  it  is  a  structural  consideration  as  the  structural  design  integration  of 
the  actuator  and  the  passive  host  structure  can  have  large  influence  on  the  energy 
efficiency  of  the  system,  but  also  it  is  a  material  consideration  as  different  actuator 
material  can  have  vastly  different  energy  efficiency.  Lastly  it  is  a  practical  issue 
concerning  the  profitability  of  the  operation  and  also  potential  regulatory  matter 

So  far  along  this  research  project,  topology  optimization  has  been  employed  to  successfully 
deal  with  mostly  the  structural  consideration,  by  formulation  objective  functions  using 
SIMP-PP  the  adaptation  requirement  and  the  stiffness  requirement  were  considered 
simultaneously  and  the  result  was  an  integrated  actuation  system  where  the  actuator  acted 
not  only  as  the  actuation  source  but  also  a  source  of  structural  integrity.  The  material  and 
practical  considerations  were  dealt  with  after  the  structural  consideration,  based  on  the 
structural  topology  given  the  selection  of  the  actuator  system  was  done  according  to  the  size 
and  control  options  available  in  a  lab  environment,  and  by  material  selection  was  conducted 
based  on  the  yield  strain  of  potential  candidates  after  FEA  analysis  of  the  topology  revealed 
the  level  of  strain  necessary  to  meet  adaptation  requirement.  This  design  behaviour,  by 
tackling  one  specific  consideration  first  and  design  the  system  based  upon  the  first 
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consideration,  then  enhance  the  design  further  by  tackling  the  remaining  two  considerations 
using  whatever  option  made  available  arising  from  the  initial  design  -  has  so  far  proven  to  be 
sufficient,  but  non-ideal  due  to  the  fact  a  large  number  of  important  parameter  are  locked 
down  after  the  initial  design,  which  limits  the  optimization  potential  of  the  final  adaptive 
system.  This  matter  should  be  a  focus  in  the  future  work  ahead  for  adaptive  airfoil  research. 
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